* Evaluate recovered beliefs from GSU risk and beliefs data on core financial literacy */

* log file
capture log close _all
log using "Evaluate Recovered Beliefs on Core Literacy Questions.log", replace name(recover)

* configure Stata
capture: version 18
set processors 4
set more off
set scheme s1color
set seed 987654321
timer clear 1
timer on 1

* packages to get
* grc1leg2

* tasks
global doIND "y"
global doPOOL "y"
global doREPORT "y"

* use reports rather than recovered beliefs (for Raven)
global useREPORTS "n"

* do risk display or not
global doRiskDisplay "n"

* number of subjects -- actually only 777-1 with matched risk preferences
local Nsub = 777

* number of questions
local Nq = 2

* name of file stem
local belief "literacy"

* graphics font
graph set window fontface "Candara"

* create directories needed
capture: mkdir figures

* tell us what version ran
about


/* belief questions for session 1

"fin1","Suppose you had $100 in a savings account and the interest rate was 2 percent per year. After 5 years, how much do you think you would have in the account if you left the money to grow?","$92","$94","$96","$98","$100","$102","$104","$106","$108","$110",25,25
"fin7","Imagine that you have $100 in a savings account and the annual interest rate on your savings account was 1 percent per year, and annual inflation was 2 percent per year. After one year, how much purchasing power would you have on the initial $100?","$95","$96","$97","$98","$99","$100","$101","$102","$103","$104",25,25

*/

* define the financial literacy questions
local ques1 = "the Savings Question"
local ques2 = "the Inflation Question"

local ques1_ = "Savings"
local ques2_ = "Inflation"

local ques3_ = "{bf:All 15 Questions}"

* the true posterior
local true1 = "$110.41"
local true2 = "$98.98"

local true1v = 110.41
local true2v = 98.98

* the true posterior as the bin number
local true1b = 10
local true2b = 5

* the true posterior as the bin number with the new labels
local true1b_new = 5
local true2b_new = 2

* location of box for title in detailed pictures
local boxloc1 = 1
local boxloc2 = 1

* location of box for legend in comparison pictures
local legloc1 = 1
local legloc2 = 1

local female0 "Men"
local female1 "Women"

* check if question was asked
qui {
use core_literacy, clear
generate int question=.
replace question=1 if qid=="fin1"
replace question=2 if qid=="fin7"
egen int sid = group(id)
forvalues s=1/`Nsub' {
forvalues q=1/`Nq' {
	count if question == `q' & sid == `s'
	local count`s'_`q' = r(N)
}
}
}

* missing risk preferences for subject sid = 182
foreach s in 182  {
forvalues q=1/`Nq' {
	local count`s'_`q' = 0
}
}

* evaluate over subjects and questions
if "$doIND" == "y" {

qui {
noi: di "Evaluating beliefs of `Nsub' subjects..."

* main calculations
forvalues s=1/`Nsub' {
forvalues q=1/`Nq' {

* number of QSR bins
local bin = 10

* skip if missing
if `count`s'_`q'' > 0 {

* get the recovered beliefs
use recovered/`belief'_s`s'_q`q', clear
summ

* label for mcmc sample
rename mcmc mcmc_size
generate long mcmc = _n
label variable mcmc "MCMC sample number"

* check beliefs sum to 1
generate check_belief = 0
forvalues x=1/`bin' {
	replace check_belief = check_belief + b_`x'
}
summ check_belief
rename check_belief normalize_b
label variable normalize_b "Normalization for recovered beliefs to sum to 1"

* get posterior predictive mean and median for each bin
generate sum_p_mean = 0
generate sum_p_median = 0
forvalues x=1/`bin' {
	summ b_`x', detail
	local m = r(mean)
	generate p_mean_`x' = `m'
	label variable p_mean_`x' "Posterior predictive mean for bin `x'"
	local m = r(p50)
	generate p_median_`x' = `m'
	label variable p_median_`x' "Posterior predictive median for bin `x'"
	replace sum_p_mean = sum_p_mean + p_mean_`x'
	replace sum_p_median = sum_p_median + p_median_`x'

	label variable r_`x' "Report for bin `x'"
	label variable b_`x' "Recovered belief for bin `x'"
}
* normalize posterior mean and median
forvalues x=1/`bin' {
	replace p_mean_`x' = p_mean_`x'/sum_p_mean
	replace p_median_`x' = p_median_`x'/sum_p_median
}
summ p_mean_*
summ p_median_*

rename sum_p_mean normalize_mean
label variable normalize_mean "Normalization for recovered posterior predictive mean to sum to 1"
rename sum_p_median normalize_median
label variable normalize_median "Normalization for recovered posterior predictive median to sum to 1"

* illustrate
des mcmc r_* b_* normalize_b normalize_mean p_mean_* normalize_median p_median_*
list mcmc r_* r eta phi in 1/5, noobs
list mcmc b_* normalize_b in 1/5, noobs
list mcmc normalize_mean p_mean_* in 1/5, noobs
list mcmc normalize_median p_median_* in 1/5, noobs

* summarize risk preference parameters and display
foreach x in r eta phi {
	summ `x'
	local `x' = r(mean)
}

* maximum likelihood parameters and beliefs: only for Stata BHM, not Stan BHM
rename _loglikelihood ll
summ ll
local maxLL = r(max)
di "Maximum LL is `maxLL'"
summ ll, detail

* get the second highest LL value, and use it to get the maxLL (avoids rounding error)
sort ll
summ ll
local nobs = r(N)
list ll if _n == `nobs'-1
list ll if _n == `nobs'
summ ll if _n == `nobs'-1
local maxLL = r(mean)
di "Maximum LL is strictly greater then `maxLL', and here it is..."
summ ll if ll > `maxLL'

* now get the ML estimates
list r eta phi ll if ll == `maxLL', noobs
foreach x in r eta phi {
	summ `x' if ll > `maxLL'
	local `x'_ml = r(mean)
	di "ML estimate of `x' is ``x'_ml'"
}
forvalues x=1/`bin' {
	 summ b_`x' if ll > `maxLL'
	 generate ml_`x' = r(mean)
}
summ ml_*

* true bin
forvalues x=1/`bin' {
	* assume old label is used
	generate int true_`x' = 0
	if `true`q'b' == `x' {
		replace true_`x' = 1
	}
	* check for new label being used
	summ new_label
	local new = r(mean)
	if `new' > 0 {
		if `true`q'b' == `x' {
			replace true_`x' = 0
		}
		if `true`q'b_new' == `x' {
			replace true_`x' = 1
		}
	}
}

* save, reduce to one observation and save the summary
save recovered/`belief'_s`s'_q`q'_save, replace
keep if _n == 1
save recovered/`belief'_s`s'_q`q'_summary, replace

* display
list, noobs
reshape long r_ b_ p_mean_ p_median_ ml_ true_ post_, i(mcmc) j(bbin)
list, noobs
rename r_ report
rename b_ belief
rename p_mean_ p_mean
rename p_median_ p_median
rename ml_ ml
rename true_ true
rename post_ post

* comparison of reports and beliefs
local angle = 0
if "`q'" == "1" {
	if `new' > 0 {
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "$102" 2 "$104" 3 "$106" 4 "$108" 5 "$110" 6 "$112" 7 "$114" 8 "$116" 9 "$118" 10 "$120") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question", size(large)) saving(figures/`belief'_s`s'_q`q', replace)
		graph bar (mean) report (mean) p_mean (mean) ml, over(bbin, relabel(1 "$102" 2 "$104" 3 "$106" 4 "$108" 5 "$110" 6 "$112" 7 "$114" 8 "$116" 9 "$118" 10 "$120") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "ML Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'_, replace)
		graph bar (mean) report (mean) p_mean (mean) true, over(bbin, relabel(1 "$102" 2 "$104" 3 "$106" 4 "$108" 5 "$110" 6 "$112" 7 "$114" 8 "$116" 9 "$118" 10 "$120") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'__, replace)
	}
	else {
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "$92" 2 "$94" 3 "$96" 4 "$98" 5 "$100" 6 "$102" 7 "$104" 8 "$106" 9 "$108" 10 "$110") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question", size(large)) saving(figures/`belief'_s`s'_q`q', replace)
		graph bar (mean) report (mean) p_mean (mean) ml, over(bbin, relabel(1 "$92" 2 "$94" 3 "$96" 4 "$98" 5 "$100" 6 "$102" 7 "$104" 8 "$106" 9 "$108" 10 "$110") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "ML Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'_, replace)
		graph bar (mean) report (mean) p_mean (mean) true, over(bbin, relabel(1 "$92" 2 "$94" 3 "$96" 4 "$98" 5 "$100" 6 "$102" 7 "$104" 8 "$106" 9 "$108" 10 "$110") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'__, replace)
	}
}
if "`q'" == "2" {
	if `new' > 0 {
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "$98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question", size(large)) saving(figures/`belief'_s`s'_q`q', replace)
		graph bar (mean) report (mean) p_mean (mean) ml, over(bbin, relabel(1 "$98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "ML Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'_, replace)
		graph bar (mean) report (mean) p_mean (mean) true, over(bbin, relabel(1 "$98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'__, replace)
	}
	else {
		graph bar (mean) report (mean) p_mean, over(bbin, relabel(1 "$95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question", size(large)) saving(figures/`belief'_s`s'_q`q', replace)
		graph bar (mean) report (mean) p_mean (mean) ml, over(bbin, relabel(1 "$95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "ML Belief") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'_, replace)
		graph bar (mean) report (mean) p_mean (mean) true, over(bbin, relabel(1 "$95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(horizontal)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Bayesian Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("Beliefs of Subject #`s' for the {it:`ques`q'_'} Question" , size(large)) saving(figures/`belief'_s`s'_q`q'__, replace)
	}
}
graph export figures/`belief'_s`s'_q`q'.png, replace

* close loop over skip for missing questions
}

* close loop over questions
}

* do risk preferences or not
if "doRiskDisplay" == "y" {

* Bayesian posterior mean risk preference parameters
clear
set obs 1

local r_ = string(`r', "%3.2f")
local phi_ = string(`phi', "%3.2f")
local eta_ = string(`eta', "%3.2f")

* initial picture of pwf
graph twoway function y = (exp( (-`eta') * (-ln(x) )^`phi')), clpattern(dash) clwidth(thick) || function y=x, clpattern(solid) legend(off) ytitle("{&omega}(p)", size(large) orientation(horizontal)) xtitle("p", size(large)) ylabel(0(.25)1, nogrid angle(horizontal)) xlabel(0(.25)1, nogrid angle(horizontal)) xline(0.5, lpattern(dash) lwidth(thin) lcolor(red*0.5)) xsize(1) ysize(1) title("r=`r_'" "{&eta}=`eta_' {&phi}=`phi_'", pos(11) ring(0) box j(center) size(medsmall)) subtitle("BHM Estimates", pos(5) ring(0) box j(center) size(medsmall)) saving(tmp1G, replace)

* now calculate the decision weights for uniform-probability lottery

* look at probability weights for various settings
foreach ncases in 7 6 5 4 3 2 {

noi {

local inc=100/`ncases'

forvalues x=0/100 {
    generate p`x'=.
    generate prob`x'=.
}

forvalues x=0(`inc')100 {
    local xx=round(`x',1)
    local xxx=`x'/100
    replace p`xx' = (exp( (-`eta') * (-ln(`xxx') )^`phi'))
    replace prob`xx'=`x'/100
}
replace p0=0
replace p100=1
replace prob0=0
replace prob100=1

summ

generate row=1
reshape long p prob, i(row) j(pr)
drop if prob==.

generate p_inc=p[_n]-p[_n-1]
drop if p_inc==.

egen prizesG=seq(), from(`ncases') to(1)
egen prizesL=seq(), from(1) to(`ncases')
generate ncases=`ncases'

}

list prizesG prizesL prob p p_inc ncases, noobs
save prob`ncases', replace

drop if prizesG>1
drop row p prob p_inc ncases prizesG prizesL pr
set obs 1

}

* collate
use prob2, clear
foreach ncases in 5 4 3 {
    append using prob`ncases'
    erase prob`ncases'.dta
}
erase prob2.dta

tab ncases
tab prizesG

summ ncases
local Nprizes=r(max)
drop if prizesG>`Nprizes'

twoway (connected p_inc prizesG if ncases==4, mcolor(black) lcolor(black)) (connected p_inc prizesG if ncases==3, mcolor(black) lcolor(black)) (connected p_inc prizesG if ncases==2, mcolor(black) lcolor(black)), ytitle("Decision" "Weight", size(large) orientation(horizontal)) ylabel(0(.25)1, angle(horizontal)) yline(0.5, lpattern(dash) lwidth(thin) lcolor(gray)) yline(0.333, lpattern(dash) lwidth(thin) lcolor(gray)) yline(0.25, lpattern(dash) lwidth(thin) lcolor(gray)) xtitle(Prize (Worst to Best), size(large)) xlabel(1(1)4) legend(off) saving(tmp2G, replace)

gr combine tmp1G.gph tmp2G.gph, title("Posterior Mean Estimates of Prelec Probability Weighting" "and Implied Decision Weights for Subject #`s'", size(large)) subtitle("Based on equi-probable reference lotteries for 2, 3 or 4 prizes", size(medsmall)) imargin(tiny) saving(p_mean_s`s', replace)


* ML parameters
clear
local r = `r_ml'
local eta = `eta_ml'
local phi = `phi_ml'

local r_ = string(`r', "%3.2f")
local phi_ = string(`phi', "%3.2f")
local eta_ = string(`eta', "%3.2f")

* initial picture of pwf
graph twoway function y = (exp( (-`eta') * (-ln(x) )^`phi')), clpattern(dash) clwidth(thick) || function y=x, clpattern(solid) legend(off) ytitle("{&omega}(p)", size(large) orientation(horizontal)) xtitle("p", size(large)) ylabel(0(.25)1, nogrid angle(horizontal)) xlabel(0(.25)1, nogrid angle(horizontal)) xline(0.5, lpattern(dash) lwidth(thin) lcolor(red*0.5)) xsize(1) ysize(1) title("r=`r_'" "{&eta}=`eta_' {&phi}=`phi_'", pos(11) ring(0) box j(center) size(medsmall)) subtitle("ML Estimates", pos(5) ring(0) box j(center) size(medsmall)) saving(tmp1G_ml, replace)

* now calculate the decision weights for uniform-probability lottery

* look at probability weights for various settings
foreach ncases in 7 6 5 4 3 2 {

qui {

local inc=100/`ncases'

forvalues x=0/100 {
    generate p`x'=.
    generate prob`x'=.
}

forvalues x=0(`inc')100 {
    local xx=round(`x',1)
    local xxx=`x'/100
    replace p`xx' = (exp( (-`eta') * (-ln(`xxx') )^`phi'))
    replace prob`xx'=`x'/100
}
replace p0=0
replace p100=1
replace prob0=0
replace prob100=1

summ
generate row=1
reshape long p prob, i(row) j(pr)
drop if prob==.

generate p_inc=p[_n]-p[_n-1]
drop if p_inc==.

egen prizesG=seq(), from(`ncases') to(1)
egen prizesL=seq(), from(1) to(`ncases')
generate ncases=`ncases'

}

list prizesG prizesL prob p p_inc ncases, noobs
save prob`ncases', replace

drop if prizesG>1
drop row p prob p_inc ncases prizesG prizesL pr
set obs 1

}

* collate
use prob2, clear
foreach ncases in 5 4 3 {
    append using prob`ncases'
    erase prob`ncases'.dta
}
erase prob2.dta

tab ncases
tab prizesG

summ ncases
local Nprizes=r(max)
drop if prizesG>`Nprizes'

twoway (connected p_inc prizesG if ncases==4, mcolor(black) lcolor(black)) (connected p_inc prizesG if ncases==3, mcolor(black) lcolor(black)) (connected p_inc prizesG if ncases==2, mcolor(black) lcolor(black)), ytitle("Decision" "Weight", size(large) orientation(horizontal)) ylabel(0(.25)1, angle(horizontal)) yline(0.5, lpattern(dash) lwidth(thin) lcolor(gray)) yline(0.333, lpattern(dash) lwidth(thin) lcolor(gray)) yline(0.25, lpattern(dash) lwidth(thin) lcolor(gray)) xtitle(Prize (Worst to Best), size(large)) xlabel(1(1)4) legend(off) saving(tmp2G_ml, replace)

gr combine tmp1G_ml.gph tmp2G_ml.gph, title("Maximum Likelihood Estimates of Prelec Probability Weighting" "and Implied Decision Weights for Subject #`s'", size(large)) subtitle("Based on equi-probable reference lotteries for 2, 3 or 4 prizes", size(medsmall)) imargin(tiny) saving(ml_s`s', replace)

gr combine tmp1G.gph tmp2G.gph tmp1G_ml.gph tmp2G_ml.gph, title("Alternative Estimates of Prelec Probability Weighting" "and Implied Decision Weights for Subject #`s'", size(large)) subtitle("Based on equi-probable reference lotteries for 2, 3 or 4 prizes", size(medsmall)) imargin(medium) saving(figures/`belief'_est_s`s', replace)
gr export figures/`belief'_est_s`s'.png, replace

* clean up
erase ml_s`s'.gph
erase p_mean_s`s'.gph
erase tmp1G.gph
erase tmp2G.gph
erase tmp1G_ml.gph
erase tmp2G_ml.gph
capture: erase prob6.dta
capture: erase prob7.dta

* close loop over doing the risk preferences
}

* close loop over subjects
}

* end of qui
}

* end of $doIND
}


* collate files
if "$doPOOL" == "y" {

* qui block
qui {

* keep track of files -- summary files
local f = 0
forvalues s=1/`Nsub' {

	* skip this one sid
	if `s' ~= 182 {

		forvalues q=1/`Nq' {

			* skip if missing
			if `count`s'_`q'' > 0 {

				* get the recovered and processed beliefs
				use recovered/`belief'_s`s'_q`q'_summary, clear

				* accumulate
				local f = `f'+1
				save tmp`f', replace

			}

		}
	}
}

* append
use tmp1, clear
forvalues x=2/`f' {
	append using tmp`x'
	erase tmp`x'.dta
}
erase tmp1.dta
save `belief'_summary, replace
noi: di "* saved the pooled summary file..."

* keep track of files -- save files
local f = 0
forvalues s=1/`Nsub' {

	* skip this one sid
	if `s' ~= 182 {

		forvalues q=1/`Nq' {

			* skip if missing
			if `count`s'_`q'' > 0 {

				* get the recovered and processed beliefs
				use recovered/`belief'_s`s'_q`q'_save, clear

				* accumulate
				local f = `f'+1
				save tmp`f', replace

			}

		}
	}
}

* append
use tmp1, clear
forvalues x=2/`f' {
	append using tmp`x'
	erase tmp`x'.dta
}
erase tmp1.dta
save `belief'_save, replace
noi: di "* saved the pooled save file with all MCMC samples..."

* end of qui
}

* review
tab sid
tab qid

* end of $doPOOL
}





* PDF document
if "$doREPORT" == "y" {

* quietly
noi {

noi: di "* Generating literacy measures for the report..."

* data for pooled displays
use `belief'_summary, clear
reshape long r_ b_ p_mean_ p_median_ ml_ true_ post_, i(sid question) j(bbin)
rename r_ report
rename b_ belief
rename p_mean_ p_mean
rename p_median_ p_median
rename ml_ ml
rename true_ true
rename post_ post

* labels for bins
generate v_lo = .
generate v_hi = . 

* assign old labels to bins
replace v_lo = .         if bbin==1 & qid=="fin1"
replace v_hi = 93        if bbin==1 & qid=="fin1"
replace v_lo = 93.01     if bbin==2 & qid=="fin1"
replace v_hi = 95        if bbin==2 & qid=="fin1"
replace v_lo = 95.01     if bbin==3 & qid=="fin1"
replace v_hi = 97        if bbin==3 & qid=="fin1"
replace v_lo = 97.01     if bbin==4 & qid=="fin1"
replace v_hi = 99        if bbin==4 & qid=="fin1"
replace v_lo = 99.01     if bbin==5 & qid=="fin1"
replace v_hi = 101       if bbin==5 & qid=="fin1"
replace v_lo = 101.01    if bbin==6 & qid=="fin1"
replace v_hi = 103       if bbin==6 & qid=="fin1"
replace v_lo = 103.01    if bbin==7 & qid=="fin1"
replace v_hi = 105       if bbin==7 & qid=="fin1"
replace v_lo = 105.01    if bbin==8 & qid=="fin1"
replace v_hi = 107       if bbin==8 & qid=="fin1"
replace v_lo = 107.01    if bbin==9 & qid=="fin1"
replace v_hi = 109       if bbin==9 & qid=="fin1"
replace v_lo = 109.01    if bbin==10 & qid=="fin1"
replace v_hi = .         if bbin==10 & qid=="fin1"

* assign new labels to bins. 
replace v_lo = .		if bbin==1 & qid=="fin1" & new_label==1
replace v_hi = 103		if bbin==1 & qid=="fin1" & new_label==1
replace v_lo = 103.01		if bbin==2 & qid=="fin1" & new_label==1
replace v_hi = 105		if bbin==2 & qid=="fin1" & new_label==1
replace v_lo = 105.01		if bbin==3 & qid=="fin1" & new_label==1
replace v_hi = 107		if bbin==3 & qid=="fin1" & new_label==1
replace v_lo = 107.01		if bbin==4 & qid=="fin1" & new_label==1
replace v_hi = 109		if bbin==4 & qid=="fin1" & new_label==1
replace v_lo = 109.01		if bbin==5 & qid=="fin1" & new_label==1
replace v_hi = 111		if bbin==5 & qid=="fin1" & new_label==1
replace v_lo = 111.01		if bbin==6 & qid=="fin1" & new_label==1
replace v_hi = 113		if bbin==6 & qid=="fin1" & new_label==1
replace v_lo = 113.01		if bbin==7 & qid=="fin1" & new_label==1
replace v_hi = 115		if bbin==7 & qid=="fin1" & new_label==1
replace v_lo = 115.01		if bbin==8 & qid=="fin1" & new_label==1
replace v_hi = 117		if bbin==8 & qid=="fin1" & new_label==1
replace v_lo = 117.01		if bbin==9 & qid=="fin1" & new_label==1
replace v_hi = 119		if bbin==9 & qid=="fin1" & new_label==1
replace v_lo = 119.01		if bbin==10 & qid=="fin1" & new_label==1
replace v_hi = .		if bbin==10 & qid=="fin1" & new_label==1

* assign old labels to bins
replace v_lo = .         if bbin==1 & qid=="fin7"
replace v_hi = 95.5      if bbin==1 & qid=="fin7"
replace v_lo = 95.51     if bbin==2 & qid=="fin7"
replace v_hi = 96.5      if bbin==2 & qid=="fin7"
replace v_lo = 96.51     if bbin==3 & qid=="fin7"
replace v_hi = 97.5      if bbin==3 & qid=="fin7"
replace v_lo = 97.51     if bbin==4 & qid=="fin7"
replace v_hi = 98.5      if bbin==4 & qid=="fin7"
replace v_lo = 98.51     if bbin==5 & qid=="fin7"
replace v_hi = 99.5      if bbin==5 & qid=="fin7"
replace v_lo = 99.51     if bbin==6 & qid=="fin7"
replace v_hi = 100.5     if bbin==6 & qid=="fin7"
replace v_lo = 100.51    if bbin==7 & qid=="fin7"
replace v_hi = 101.5     if bbin==7 & qid=="fin7"
replace v_lo = 101.51    if bbin==8 & qid=="fin7"
replace v_hi = 102.5     if bbin==8 & qid=="fin7"
replace v_lo = 102.51    if bbin==9 & qid=="fin7"
replace v_hi = 103.5     if bbin==9 & qid=="fin7"
replace v_lo = 103.51    if bbin==10 & qid=="fin7"
replace v_hi = .         if bbin==10 & qid=="fin7"

* assign new labels to bins. 
replace v_lo = .	if bbin==1 & qid=="fin7" & new_label==1
replace v_hi = 98.5	if bbin==1 & qid=="fin7" & new_label==1
replace v_lo = 98.51	if bbin==2 & qid=="fin7" & new_label==1
replace v_hi = 99.5	if bbin==2 & qid=="fin7" & new_label==1
replace v_lo = 99.51	if bbin==3 & qid=="fin7" & new_label==1
replace v_hi = 100.5	if bbin==3 & qid=="fin7" & new_label==1
replace v_lo = 100.51	if bbin==4 & qid=="fin7" & new_label==1
replace v_hi = 101.5	if bbin==4 & qid=="fin7" & new_label==1
replace v_lo = 101.51	if bbin==5 & qid=="fin7" & new_label==1
replace v_hi = 102.5	if bbin==5 & qid=="fin7" & new_label==1
replace v_lo = 102.51	if bbin==6 & qid=="fin7" & new_label==1
replace v_hi = 103.5	if bbin==6 & qid=="fin7" & new_label==1
replace v_lo = 103.51	if bbin==7 & qid=="fin7" & new_label==1
replace v_hi = 104.5	if bbin==7 & qid=="fin7" & new_label==1
replace v_lo = 104.51	if bbin==8 & qid=="fin7" & new_label==1
replace v_hi = 105.5	if bbin==8 & qid=="fin7" & new_label==1
replace v_lo = 105.51	if bbin==9 & qid=="fin7" & new_label==1
replace v_hi = 106.5	if bbin==9 & qid=="fin7" & new_label==1
replace v_lo = 106.51	if bbin==10 & qid=="fin7" & new_label==1
replace v_hi = .	if bbin==10 & qid=="fin7" & new_label==1

* get report as integer summing to 100
generate int choiceI = round(belief*100,1)
summ belief choiceI

* replace with report
if "$useREPORTS" == "y" {
	replace choiceI = round(report*100,1)
	summ report choiceI
}

* get mean and standard deviation
intreg v_lo v_hi [fweight = choiceI] if qid=="fin7", cluster(sid)
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local sd_ = string(`sd', "%4.2f")
di "Mean is `mean_', sigma is `sigma' and standard deviation is `sd_'"

* all subjects
local angle = 45

local q = 1
graph bar (mean) report (mean) p_mean (mean) true if question==`q' & new_label==0, over(bbin, relabel(1 "$102" 2 "$104" 3 "$106" 4 "$108" 5 "$110" 6 "$112" 7 "$114" 8 "$116" 9 "$118" 10 "$120") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp1, replace)
graph bar (mean) report (mean) p_mean (mean) true if question==`q' & new_label==1, over(bbin, relabel(1 "$92" 2 "$94" 3 "$96" 4 "$98" 5 "$100" 6 "$102" 7 "$104" 8 "$106" 9 "$108" 10 "$110") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp2, replace)
local q = 2
graph bar (mean) report (mean) p_mean (mean) true if question==`q' & new_label==0, over(bbin, relabel(1 "$98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp3, replace)
graph bar (mean) report (mean) p_mean (mean) true if question==`q' & new_label==1, over(bbin, relabel(1 "$95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(black) lcolor(black)) bar(2, fcolor(ltblue) lcolor(ltblue)) bar(3, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Report" 2 "Belief" 3 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(1)) saving(tmp4, replace)
grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, rows(2) imargin(small) title("Reports, Beliefs and True Answers for all Subjects", size(large) margin(small)) saving(figures/pooled, replace)
gr export figures/pooled.png, replace

* men
local f = 0

local q = 1
graph bar (mean) p_mean (mean) true if question==`q' & new_label==0 & female==`f', over(bbin, relabel(1 "$102" 2 "$104" 3 "$106" 4 "$108" 5 "$110" 6 "$112" 7 "$114" 8 "$116" 9 "$118" 10 "$120") gap(*0.5) label(angle(`angle') labsize(medsmall)))  bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp1, replace)
graph bar (mean) p_mean (mean) true if question==`q' & new_label==1 & female==`f', over(bbin, relabel(1 "$92" 2 "$94" 3 "$96" 4 "$98" 5 "$100" 6 "$102" 7 "$104" 8 "$106" 9 "$108" 10 "$110") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp2, replace)
local q = 2
graph bar (mean) p_mean (mean) true if question==`q' & new_label==0 & female==`f', over(bbin, relabel(1 "$98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp3, replace)
graph bar (mean) p_mean (mean) true if question==`q' & new_label==1 & female==`f', over(bbin, relabel(1 "$95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(1)) saving(tmp4, replace)
grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, rows(2) imargin(small) title("Beliefs and True Answers for Men", size(large) margin(small)) saving(figures/pooled_male, replace)
gr export figures/pooled_male.png, replace

* women
local f = 1
local q = 1
graph bar p_mean (mean) true if question==`q' & new_label==0 & female==`f', over(bbin, relabel(1 "$102" 2 "$104" 3 "$106" 4 "$108" 5 "$110" 6 "$112" 7 "$114" 8 "$116" 9 "$118" 10 "$120") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp1, replace)
graph bar (mean) p_mean (mean) true if question==`q' & new_label==1 & female==`f', over(bbin, relabel(1 "$92" 2 "$94" 3 "$96" 4 "$98" 5 "$100" 6 "$102" 7 "$104" 8 "$106" 9 "$108" 10 "$110") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp2, replace)
local q = 2
graph bar (mean) p_mean (mean) true if question==`q' & new_label==0 & female==`f', over(bbin, relabel(1 "$98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(11)) saving(tmp3, replace)
graph bar (mean) p_mean (mean) true if question==`q' & new_label==1 & female==`f', over(bbin, relabel(1 "$95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "Belief" 2 "Correct") rows(1) size(small) position(12) ring(1)) title("`ques`q'_'", box ring(0) pos(1)) saving(tmp4, replace)
grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, rows(2) imargin(small) title("Beliefs and True Answers for Women", size(large) margin(small)) saving(figures/pooled_female, replace)
gr export figures/pooled_female.png, replace


if "$useREPORTS" == "y" {
	local belief_text "Reported Beliefs"
}
else {
	local belief_text "Recovered Beliefs"
}

local f = 0
local q = 2
intreg v_lo v_hi [fweight = choiceI] if question==`q' & new_label==0 & female==`f', cluster(sid)
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
local sd_ = string(`sd', "%4.2f")
local ratio = (`bias'/`sd')*100
local ratio_ = string(`ratio', "%2.0f")
graph bar (mean) p_mean (mean) true if question==`q' & new_label==0 & female==`f', over(bbin, relabel(1 "< $95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("`female`f''", box ring(0) pos(11)) subtitle("True=$`true_', {&mu}=$`mean_', so |bias|=$`bias_'," " which is just `ratio_'% of {&sigma}=$`sd_'", margin(medsmall) size(medsmall)) saving(tmp5, replace)

intreg v_lo v_hi [fweight = choiceI] if question==`q' & new_label==1 & female==`f', cluster(sid)
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
local sd_ = string(`sd', "%4.2f")
local ratio = (`bias'/`sd')*100
local ratio_ = string(`ratio', "%2.0f")
graph bar (mean) p_mean (mean) true if question==`q' & new_label==1 & female==`f', over(bbin, relabel(1 "< $98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("`female`f''", box ring(0) pos(1)) subtitle("True=$`true_', {&mu}=$`mean_', so |bias|=$`bias_'," " which is just `ratio_'% of {&sigma}=$`sd_'", margin(medsmall) size(medsmall)) saving(tmp6, replace)

local f = 1
intreg v_lo v_hi [fweight = choiceI] if question==`q' & new_label==0 & female==`f', cluster(sid)
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
local sd_ = string(`sd', "%4.2f")
local ratio = (`bias'/`sd')*100
local ratio_ = string(`ratio', "%2.0f")
graph bar (mean) p_mean (mean) true if question==`q' & new_label==0 & female==`f', over(bbin, relabel(1 "< $95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("`female`f''", box ring(0) pos(11)) subtitle("True=$`true_', {&mu}=$`mean_', so |bias|=$`bias_'," " which is just `ratio_'% of {&sigma}=$`sd_'", margin(medsmall) size(medsmall)) saving(tmp7, replace)

intreg v_lo v_hi [fweight = choiceI] if question==`q' & new_label==1 & female==`f', cluster(sid)
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
local sd_ = string(`sd', "%4.2f")
local ratio = (`bias'/`sd')*100
local ratio_ = string(`ratio', "%2.0f")
graph bar (mean) p_mean (mean) true if question==`q' & new_label==1 & female==`f', over(bbin, relabel(1 "< $98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("`female`f''", box ring(0) pos(1)) subtitle("True=$`true_', {&mu}=$`mean_', so |bias|=$`bias_'," " which is just `ratio_'% of {&sigma}=$`sd_'", margin(medsmall) size(medsmall)) saving(tmp8, replace)

grc1leg2 tmp5.gph tmp6.gph tmp7.gph tmp8.gph, cols(2) imargin(small) title("Bias and Confidence in the Literacy of" "Men and Women on `ques`q''", size(large) margin(medsmall)) saving(figures/pooled_inflation, replace)

local fig = 17
if "$useREPORTS" == "n" {
	grc1leg2 tmp5.gph tmp6.gph tmp7.gph tmp8.gph, cols(2) imargin(small) title("Figure `fig': Bias and Confidence in the Literacy of" "Men and Women on `ques`q''", size(large) margin(medsmall)) saving(raven_`fig'.gph, asis replace)
	gr export figures/pooled_inflation.png, replace
	grc1leg2 tmp5.gph tmp6.gph tmp7.gph tmp8.gph, cols(2) imargin(small) title("Bias and Confidence in the Literacy of" "Men and Women on `ques`q''", size(large) margin(medsmall)) saving(raven_`fig'_ppt.gph, asis replace)
}
else {
	grc1leg2 tmp5.gph tmp6.gph tmp7.gph tmp8.gph, cols(2) imargin(small) title("Figure B`fig': Bias and Confidence in the Literacy of" "Men and Women on `ques`q''", size(large) margin(medsmall)) saving(raven_B`fig'.gph, asis replace)
	grc1leg2 tmp5.gph tmp6.gph tmp7.gph tmp8.gph, cols(2) imargin(small) title("Bias and Confidence in the Literacy of" "Men and Women on `ques`q''", size(large) margin(medsmall)) saving(raven_B`fig'_ppt.gph, asis replace)
}

* get individual pictures to illustrate
local q = 2

* perfectly iliterate
local s = 756 
intreg v_lo v_hi [fweight = choiceI] if question==`q' & sid==`s'
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
local sd_ = string(`sd', "%4.2f")
local ratio = (`bias'/`sd')*100
local ratio_ = string(`ratio', "%2.0f")
graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("|bias|=$`bias_', with no confidence at all", margin(medsmall) size(medsmall)) saving(tmp3, replace)

* some bias
local s = 11 
intreg v_lo v_hi [fweight = choiceI] if question==`q' & sid==`s'
nlcom (sigma: exp([lnsigma]_cons))
local sigma = exp(e(b)[1,2])
local sd = `sigma'
local mean = e(b)[1,1]
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
local sd_ = string(`sd', "%4.2f")
local ratio = (`bias'/`sd')*100
local ratio_ = string(`ratio', "%2.0f")
*graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("|bias|=$`bias_', `ratio_'% of {&sigma}", margin(medsmall) size(medsmall)) saving(tmp2, replace)
graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("|bias|=$`bias_', `ratio_'% of {&sigma}", margin(medsmall) size(medsmall)) saving(tmp2, replace)

* dogmatic bias (intreg bombs)
local s = 12
local mean = 101
local mean_ = string(`mean', "%5.2f")
local true_ = string(`true`q'v', "%5.2f")
local bias = abs(`true`q'v' - `mean')
local bias_ = string(`bias', "%3.2f")
*graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("|bias|=$`bias_', with complete confidence", margin(medsmall) size(medsmall)) saving(tmp4, replace)
graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("|bias|=$`bias_', with complete confidence", margin(medsmall) size(medsmall)) saving(tmp4, replace)

* perfectly literate 
local s = 10
*graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $98" 2 "$99" 3 "$100" 4 "$101" 5 "$102" 6 "$103" 7 "$104" 8 "$105" 9 "$106" 10 "$107 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("Zero bias, with complete confidence", margin(medsmall) size(medsmall)) saving(tmp1, replace)
graph bar (mean) p_mean (mean) true if question==`q' & sid==`s', over(bbin, relabel(1 "< $95" 2 "$96" 3 "$97" 4 "$98" 5 "$99" 6 "$100" 7 "$101" 8 "$102" 9 "$103" 10 "$104 +") gap(*0.5) label(angle(`angle') labsize(medsmall))) bar(1, fcolor(ltblue) lcolor(ltblue)) bar(2, fcolor(red*0.5) lcolor(red*0.5)) bargap(0) outergap(0) ytitle("Probability", size(medlarge) margin(small) orientation(vertical)) ylabel(0(0.25)1, labsize(medsmall) angle(horizontal)) legend(order(1 "`belief_text'" 2 "Correct") rows(1) size(small) position(6) ring(1)) title("Subject #`s'", box ring(0) pos(1) size(medium)) subtitle("Zero bias, with complete confidence", margin(medsmall) size(medsmall)) saving(tmp1, replace)

grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, cols(2) imargin(small) title("Literacy and Confidence: The Good, the Bad, and the Ugly", size(large) margin(medsmall)) saving(figures/examples_inflation, replace)

local fig = 18
if "$useREPORTS" == "n" {
	grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, cols(2) imargin(small) title("Figure `fig': Literacy and Confidence:" "The Good, the Bad, and the Ugly", size(large) margin(medsmall)) saving(raven_`fig'.gph, asis replace)
	gr export figures/examples_inflation.png, replace
	grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, cols(2) imargin(small) title("Literacy and Confidence:" "The Good, the Bad, and the Ugly", size(large) margin(medsmall)) saving(raven_`fig'_ppt.gph, asis replace)
}
else {
	grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, cols(2) imargin(small) title("Figure B`fig': Literacy and Confidence:" "The Good, the Bad, and the Ugly", size(large) margin(medsmall)) saving(raven_B`fig'.gph, asis replace)
	grc1leg2 tmp1.gph tmp2.gph tmp3.gph tmp4.gph, cols(2) imargin(small) title("Literacy and Confidence:" "The Good, the Bad, and the Ugly", size(large) margin(medsmall)) saving(raven_B`fig'_ppt.gph, asis replace)
}

forvalues x=1/8 {
	capture: erase tmp`x'.gph
}

* the data for the hypothesis tests
use `belief'_save, clear

* literacy using beliefs
generate literate = 0
forvalues x=1/10 {
	replace literate = b_`x' if true_`x'==1
}

summ literate if question==1
local x = r(mean)*100
local lit1 = string(`x', "%2.0f")

summ literate if question==2
local x = r(mean)*100
local lit2 = string(`x', "%2.0f")

summ literate if question==1 & female==1
local x = r(mean)*100
local lit1f = string(`x', "%2.0f")
summ literate if question==2 & female==1
local x = r(mean)*100
local lit2f = string(`x', "%2.0f")

summ literate if question==1 & female==0
local x = r(mean)*100
local lit1m = string(`x', "%2.0f")
summ literate if question==2 & female==0
local x = r(mean)*100
local lit2m = string(`x', "%2.0f")


* literacy using reports
replace literate = 0
forvalues x=1/10 {
	replace literate = r_`x' if true_`x'==1 & mcmc==1
}

summ literate if question==1 & mcmc==1
local x = r(mean)*100
local rlit1 = string(`x', "%2.0f")

summ literate if question==2 & mcmc==1
local x = r(mean)*100
local rlit2 = string(`x', "%2.0f")

summ literate if question==1 & female==1 & mcmc==1
local x = r(mean)*100
local rlit1f = string(`x', "%2.0f")
summ literate if question==2 & female==1 & mcmc==1
local x = r(mean)*100
local rlit2f = string(`x', "%2.0f")

summ literate if question==1 & female==0 & mcmc==1
local x = r(mean)*100
local rlit1m = string(`x', "%2.0f")
summ literate if question==2 & female==0 & mcmc==1
local x = r(mean)*100
local rlit2m = string(`x', "%2.0f")

* percent of new labels
summ new_label if question==1 & mcmc==1
local x = 100-(r(mean)*100)
local old1 = string(`x', "%2.0f")
local x = r(mean)*100
local new1 = string(`x', "%2.0f")

summ new_label if question==2 & mcmc==1
local x = 100-(r(mean)*100)
local old2 = string(`x', "%2.0f")
local x = r(mean)*100
local new2 = string(`x', "%2.0f")

* number of subjects in each question
tab sid if question==1 & mcmc==1
local x = r(r)
local Nsub1 = string(`x', "%3.0f")
tab sid if question==2 & mcmc==1
local x = r(r)
local Nsub2 = string(`x', "%3.0f")

* adjust Nsub for 1 MIA
local Nsub_ = `Nsub' - 1

* display a subset of all possible subjects
local Nsub_report = `Nsub'
local Nsub_report = 100

* start a PDF
	noi: di "* Generating PDF report..."

	putpdf begin, font("Candara",12) margin(left,1) margin(right,1)

	putpdf paragraph, font("",14) halign(center)
	if "$useREPORTS" == "y" {
		putpdf text ("Reported Beliefs and Core Financial Literacy"), bold linebreak(2)
	}
	else {
		putpdf text ("Recovered Beliefs and Core Financial Literacy"), bold linebreak(2)
	}

	local d = c(current_date)
	local t = c(current_time)
	putpdf paragraph, font("",12) halign(center)
	putpdf text ("`d'   `t'"), linebreak(2)

	putpdf paragraph, font("",12) halign(left)

	putpdf text ("        There have been `Nsub_' subjects, spanning sessions in 2013, 2014 and 2016, each reporting their beliefs about 2 core financial literacy questions. ")
	putpdf text ("One question, asked of `Nsub1' subjects, is: ")
	putpdf text ("Suppose you had $100 in a savings account and the interest rate was 2 percent per year. After 5 years, how much do you think you would have in the account if you left the money to grow? "), italic
	putpdf text ("The correct answer is `true1'. The other question, asked of `Nsub2' subjects, is: ")
	putpdf text ("Imagine that you have $100 in a savings account and the annual interest rate on your savings account was 1 percent per year, and annual inflation was 2 percent per year. After one year, how much purchasing power would you have on the initial $100? "), italic
	putpdf text ("The correct answer is `true2'. The subjects that answered the first question also answered the second question. These questions are based on two of the core financial literacy questions posed by Lusardi and Mitchell. They asked the questions in ")
	putpdf text ("a multiple-choice format, with no incentives for the correct answer. ")

	putpdf text (" "), linebreak(2)
	putpdf text ("        The belief interface we used employs a Quadratic Scoring Rule (QSR) ")
	putpdf text ("with parameters chosen to allow a perfect allocation of 100 tokens to the correct interval to be rewarded with $50 if that questions was chosen for payment. ")
	putpdf text ("The QSR reports were defined over 10 bin intervals, and one interval contained the correct answer. ")
	putpdf text ("The risk preferences were estimated with a Bayesian Hierarchical Model (BHM), assuming a RDU model of risk preferences with Prelec probability weighting. ")
	if "$useREPORTS" == "y" {
		putpdf text ("Beliefs are assumed to be the same as reports in these displays. ")
	}
	else {
		putpdf text ("Beliefs were recovered from these reports using the individual risk estimates from the BHM for each subject, and assuming that they responded rationally when reporting beliefs conditional on those risk preferences. ")
	}

	putpdf text (" "), linebreak(2)
	putpdf text ("        A treatment was to shift the labels so that the correct answer was in a different belief report bin interval. ")
	putpdf text ("For the first question the correct bin was #10 for `old1'% of the subjects, and #5 for `new1'% of the subjects. ")
	putpdf text ("For the second question the correct bin was #5 for `old2'% of the subjects, and #2 for `new2'% of the subjects. Results on literacy are reported across this treatment. ")

	putpdf text (" "), linebreak(2)
	putpdf text ("        For ")
	putpdf text ("`ques1'"), bold
	putpdf text (", the overall literacy in terms of reports was `rlit1'%. The literacy for men was `rlit1m'% and the literacy for women was `rlit1f'%  ")
	putpdf text ("For `ques1', the overall literacy in terms of beliefs was `lit1'%. The literacy for men was `lit1m'% and the literacy for women was `lit1f'%  ")

	putpdf text (" "), linebreak(2)
	putpdf text ("        For ")
	putpdf text ("`ques2'"), bold
	putpdf text (", the overall literacy in terms of reports was `rlit2'%. The literacy for men was `rlit2m'% and the literacy for women was `rlit2f'%  ")
	putpdf text ("For `ques2', the overall literacy in terms of beliefs was `lit2'%. The literacy for men was `lit2m'% and the literacy for women was `lit2f'%  ")

	putpdf text (" "), linebreak(2)
	putpdf text ("        The first two displays are the ones intended for the revision of the Raven paper, focussing solely on the second question (the first question gets uniformly high literacy). ")
	putpdf text ("The displays then documents pooled reports, beliefs and true answers for all subjects, then pooled beliefs and true answers for men, and pooled beliefs and true answers for women. It then documents reports, beliefs, and true answers for each individual subject for the second question. ")
	if `Nsub_report' < `Nsub' {
		putpdf text ("For purposes of illustration we only report individual results for `Nsub_report' subjects, and then only if they answered both questions. ")
	}

	* list the last page of text
	local page = 2

	* pooled beliefs
	putpdf pagebreak
	putpdf paragraph, halign(right) font("Candara",10)
	local page = `page' + 1
	putpdf text ("Page `page'"), linebreak(1)
	putpdf paragraph, halign(center)
	putpdf image figures/pooled_inflation.png, linebreak width(6)
	putpdf image figures/examples_inflation.png, linebreak width(6)

	putpdf pagebreak
	putpdf paragraph, halign(right) font("Candara",10)
	local page = `page' + 1
	putpdf text ("Page `page'"), linebreak(1)
	putpdf paragraph, halign(center)
	putpdf image figures/pooled.png, linebreak width(6)

	putpdf pagebreak
	putpdf paragraph, halign(right) font("Candara",10)
	local page = `page' + 1
	putpdf text ("Page `page'"), linebreak(1)
	putpdf paragraph, halign(center)
	putpdf image figures/pooled_male.png, linebreak width(6)
	putpdf image figures/pooled_female.png, linebreak width(6)

	* individual beliefs
	use `belief'_summary, clear
	forvalues s = 1/`Nsub_report' {
		* skip this one sid
		if `s' ~= 182 {
			tab question if sid==`s' & question==2
			local n = r(r)
			* only if question 2 asked
			if `n'>0 {
				putpdf pagebreak
				putpdf paragraph, halign(right) font("Candara",10)
				local page = `page' + 1
				putpdf text ("Page `page'"), linebreak(1)
				putpdf paragraph, halign(center)
				putpdf image figures/`belief'_s`s'_q2.png, linebreak width(6)
			}
		}
	}

* save the PDF report
if "$useREPORTS" == "y" {
	putpdf save "Core Financial Literacy Measures from Reported Beliefs.pdf", replace
}
else {
	putpdf save "Core Financial Literacy Measures from Recovered Beliefs.pdf", replace
}

* end of quietly
}

* end of $doREPORT
}

* time taken overall
timer off 1
timer list
local secs = r(t1)
local mins = `secs'/60
local hrs = `mins'/60
local secs_ = string(`secs', "%10.0f")
local mins_ = string(`mins', "%4.1f")
local hrs_ = string(`hrs', "%4.2f")
di "Complete calculations took `secs_' seconds, `mins_' minutes, or `hrs_' hours."

log close recover
